knitr::opts_chunk$set(
    echo = TRUE,
    warning = FALSE,
    comment = "##",
    prompt = TRUE,
    tidy = TRUE,
    tidy.opts = list(width.cutoff = 75),
    fig.path = "img/"
)

Before we take you on our expedition to learn how to analyze home range data, be sure to install the following packages in R: {maptools}, {sp}, {rgdal}, {zoom}, {adehabitatHR}.

> library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
> library(sp)
> library(rgdal)
## rgdal: version: 1.2-4, (SVN revision 643)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
##  Path to GDAL shared files: C:/Users/Silvy/Documents/R/win-library/3.3/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: C:/Users/Silvy/Documents/R/win-library/3.3/rgdal/proj
##  Linking to sp version: 1.2-3

We will be looking at data from three groups if titi monkeys (Callicebus discolor) that are studied at the Tiputini Biodiversity Station in Ecuador. Let’s first load in a polygon of the trail sytem at research station. This can be done with the readShapeLines() function in the maptools package.

> tbs <- readShapeLines("C:/Users/Silvy/Documents/Austin/Tiputini/GIS/trails.shp")
> plot(tbs)

One of the complex parts of working with these file types is dealing with coordinate reference systems. This is largely due to the fact that Geographical Coordinate Systems and Projected Coordinate Systems are different in how they project the bumpy ellipsoid that is our planet onto a plane.

So once you have loaded in a Shapefile, it is important to first find out what coordinate system is used, if any. This is done easily with with proj4string commands and functions in the package sp.

> tbs@proj4string  # This command shows us that there is no Coordinate System assigned to the TBS shapefile. 
## CRS arguments: NA

To be sure there are no coordinate data related to the file in any way, run the following code. Even though it seems there is no chosen projection, in this case the new projection will not be used correctly without using this command.

> proj4string(tbs) <- NA_character_  # remove Coordinate Reference System information from TBS shapefile.

Now we can specify our new coordinate system, which will be the Projected Coordinate System UTM (Universal Transverse Mercator). This projects the earth on a flat surface and then uses x and y coordinates to locate points, just like in the figure below:

> proj4string(tbs) <- CRS("+proj=utm +zone=18 +south +datum=WGS84")  #Setting UTM WGS84 as our Coortinate Reference System (CRS).
> plot(tbs)  #Tadaaa!

Now that we have defined the coordinate system of the shapefile, we can pull in the ranging data collected from the titi monkeys between July 2014 to June 2015. These data have been cleaned in advance and stored in a .csv file on GitHub.

> library(curl)
> f <- curl("https://raw.githubusercontent.com/Callicebus/vignette/master/GPScoordinates.csv")
> points <- read.csv(f)
> head(points)
##   Obs.Sample.ID     Date Month        Group     UTMX    UTMY
## 1       OS92719 4-Jul-14     7 Callicebus L 371667.0 9929338
## 2       OS92719 4-Jul-14     7 Callicebus L 371647.7 9929378
## 3       OS92719 4-Jul-14     7 Callicebus B 371664.1 9929523
## 4       OS92719 4-Jul-14     7 Callicebus B 371747.5 9929529
## 5       OS92719 4-Jul-14     7 Callicebus B 371711.8 9929542
## 6       OS92719 4-Jul-14     7 Callicebus B 371700.1 9929542

As you can see, the coordinates in this file are already in the UTM projection, but we will have to tell R that. This is also done through commands in the sp package. First, you’ll have to convert the .csv file to a dataframe for spatial points, and then set the right coordinate system.

> spdf <- SpatialPointsDataFrame(coords = c(points[5], points[6]), data = points[4])  #Column 5 has the X-coordinates, column 6 has the Y-coordinates, and column 4 has important attribute data.
> str(spdf)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  1719 obs. of  1 variable:
##   .. ..$ Group: Factor w/ 3 levels "Callicebus B",..: 3 3 1 1 1 1 1 1 2 2 ...
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:1719, 1:2] 371667 371648 371664 371747 371712 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "UTMX" "UTMY"
##   ..@ bbox       : num [1:2, 1:2] 371376 9929241 372130 9929694
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "UTMX" "UTMY"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA
> plot(spdf, pch = 1, col = spdf$Group)

> utm <- SpatialPoints(spdf, proj4string = CRS("+proj=utm +zone=18 +south +datum=WGS84"))
> utm <- SpatialPointsDataFrame(utm, data = points[4])
> plot(tbs)
> plot(utm, pch = 20, col = utm$Group, add = T)  #Simply adding add = TRUE to your code will overlay the second plot on the first.

Now that we have all points plotted, we can draw polygons around the points of each group, representing the home ranges of the monkeys during the one-year period. We’ll do this with the help of a package build especially for home range measurements. NOTE: the functions in this package work best if you use the UTM projections!

> library(adehabitatHR)
## Loading required package: deldir
## deldir 0.1-12
## Loading required package: ade4
## Loading required package: adehabitatMA
## Loading required package: adehabitatLT
## Loading required package: CircStats
## Loading required package: MASS
## Loading required package: boot
> library(ggplot2)
> kernel <- kernelUD(utm[, 1], h = "href")
> TitiRange <- getverticeshr(kernel, 95)
> plot(tbs)
> plot(TitiRange, col = TitiRange@data$id, add = TRUE)

> F <- fortify(TitiRange)
## Regions defined for each Polygons

We can now use these polygons to estimate home range size quite simply:

> as.data.frame(TitiRange)
##                        id     area
## Callicebus B Callicebus B 6.195812
## Callicebus K Callicebus K 4.276863
## Callicebus L Callicebus L 6.643743

Let’s look only at the home range of Callicebus group L. Using our points, we can also find out which parts of their home range re used more often, and which parts are used less frequently. First, let’s subset the data so we only have the ranging points of Callicebus L.

> CallicebusL <- utm[utm@data$Group == "Callicebus L", ]
> plot(tbs)
> plot(CallicebusL, pch = 2, col = "blue", add = T)

With the use of ggplot, we can create a heat map. This map will show the density of the ranging points throughout the home range.

library(ggplot2)
Lpoints <- points[points$Group == "CallicebusL"]
hothothot <- ggplot(Lpoints, aes(x=UTMX, y=UTMY)) + geom_point() + stat_density2d(aes(fill=..density..), geom = "tile", contour = FALSE) + scale_fill_gradient2(low = "white", high = "red")
hothothot